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ABSTRACT 

The rare and conspicuous flux density variations of some radio sources (extragalactic and pulsars) for periods of weeks to months 
have been denoted Extreme Scattering Events (ESE's) by Fiedler et al. (1987). Presently, there is no astrophysical mechanism that 
satisfactorily produces this phenomenon. In this paper, we conjecture that inhomogeneities of the electronic density in the turbulent 
interstellar medium might be the origin of this phenomenon. We have tested this conjecture by a simulation of the scintillation of the 
pulsar B1937+21 at 1.4 GHz and 1.7 GHz for a period of six months. To this end, we have constructed a large square Kolmogorov 
phase screen made of 131k x 131k pixels with electron inhomogeneity scales ranging from 6 10 6 m to 10 12 m and used the Kirchhoff- 
Fresnel integral to simulate dynamic spectra of a pulsar within the framework of Physical Optics. 

The simulated light curves exhibit a 10 day long variation simultaneously at 1.41 and 1.7 GHz that is alike the "ESE" observed with 
the Nancay radiotelescope toward the pulsar B 1937+21 in October 1989. Consequently, we conclude that "ESE" toward pulsars can 
be caused naturally by the turbulence in the ionized interstellar medium instead of invoking the crossing of discrete over pressured 
ionized clouds on the line of sight as in the model of Fiedler et al. (1987). We suggest that longer events could occur in a simulation 
of scintillation, if larger electron inhomogeneities (> 10 12 m) were included in the construction of the Kolmogorov phase screen. This 
next step requires a supercomputer. 
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^ ' 1. Introduction The mechanism usually invoked for this phenomenon is scat- 
v> , tering by a plasma lens that occurs when a discrete ionized cloud 
£j . The rare and conspicuous flux density variations of some ra- crosse s the line of sight. Fiedler et al. (1987) proposed that the 
Ctf ■ dio sources for periods of weeks to months have been denoted me dium inside the cloud is highly turbulent, causing extreme 
Extreme Scattering Events (ESE's). Usually, the flux density sca ttering responsible for the flux density variation and source 
drops by at least 50 % and increases somewhat at ingress and broadening. Differently, Romani et al. (1987) proposed that the 
egress although there is a variety of possible shapes. The first me dium inside the cloud acts as a purely refractive lens pro- 
event was recognized in the direction of the quasar 0954+658 ducing extreme refraction with caustics responsible for the flux 
at 2.7 and 8.4 GHz (Fiedler et al. 1987). The first event associ- variation and source displacement during the event. The obser- 
ated with a pulsar was detected in the direction of the millisec- vations of ESE ' S in directions of quasars and pulsars and these 
ond pulsar B 1937+21 at 1.41 GHz (Cognard et al. 1993). The mo dels provide an estimation of the cloud size of 1-50 AU and 
most recent census indicates that 15 ESE's have been identified internal electronic density of 100 - 1000 e" cm 3 (Fiedler et 
in the radio light curves of 12 quasars in the entire Green Bank al 1994) Therefore, such a structure is over pressurized rel- 
Interferometer monitoring of 149 radio sources at GHz frequen- ative to the ambient medium making its lifetime as short as a 
cies between 1979 and 1996 (Lazio et al. 2001). Five ESE's have few year s. In addition, the observed rate of occurrence implies 
been identified in the direction of the pulsar B 1937 +21 between a cloud space density as high as ~ 10 6 clouds pc 3 (LRC98). 
1989 and 1996 (Lestrade, Rickett & Cognard 1998, LRC98 here- There i s n0 known mechanism that can produce so many dis- 
after) and a very long event was observed in the direction of the crete i on ized clouds in the ISM. In order to avoid this difficulty, 
pulsar J 1643- 1224 (Maitia, Lestrade & Cognard 2003). Hill et we conjecture in this paper that the Kolmogorov turbulence of 
al. (2005) also report criss-cross pattern in the dynamic spectra the ionized ISM is the natural cause for the ESE's and their oc- 
of PSR B0834+06 at 327MHz they interpret as resulting from currence is statistical 
compact scattering objects essentially stationary in the screen in 

20 years. j n order to support this proposition, we have simulated in- 
terstellar scintillation at 1.41 and 1.7 GHz toward the pulsar 
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radiotelescopeQ. The flux density of this pulsar measured at the 
Nancay radiotelescope is an average over an integration time of 
70 minutes and an integration bandwidth of 7.5 MHz centered 
around 1.41 GHz and 12.5 MHz around 1.7 GHz. These intervals 
are large enough to average several diffractive interstellar scintil- 
lation patterns in the time-frequency domain, and so essentially 
sample the refractive regime. However, some significant influ- 
ence of diffractive ISS remains. Therefore, we conducted a full 
scintillation calculation including the diffractive and refractive 
scales. The refractive time scales measured at Nangay in the di- 
rection of B1937+21 are 16+1 days and 8+1 days, at 1.41 GHz 
and 1.7 GHz respectively. Also, the modulation indexes are m r - 
0.30+0.02 and 0.50+0.04 at 1.41 GHz and 1.7 GHz respectively. 
The zero-lag correlation of the flux density between these two 
frequencies is 0.93+0.05. Detailed descriptions of these mea- 
surements are given in LRC98. 

For the simulation presented in this paper, we have used 
the Kolmogorov 3-D power spectrum of the electronic den- 
sity 5n e (x,y,z) in the interstellar medium as established locally 
(< lkpc) by Armstrong, Rickett & Spangler (1995). We have 
used the Kirchhoff-Fresnel integral to compute pulsar dynamic 
spectra in the framework of Physical Optics. In section 2, we 
present the formulae for this simulation carried out in the thin 
screen approximation, the construction of the Kolmogorov phase 
screen and the computation of the diffraction pattern (see also 
Hamidouche 2003). In section 3, we construct a large square 
Kolmogorov phase screen of 131kx 131k pixels. In section 4, 
we make a series of tests to empirically determine how to set 
the parameters of the Kirchoff-Fresnel integral. In section 5, we 
present the result of the simulation in the direction of B 1937+21 
simulated over a period of 6 months. 

2. Scattering of a Point-Like Source by the Ionized 
Interstellar Medium 

The historical problem of diffraction have been cast in mod- 
ern mathematical terms by Born and Wolf (1999) and Goodman 
(1968). A particularly clear presentation of this problem applied 
to interstellar optics in the radio domain is in Narayan (1992) 
and Gwinn et al. (1998). The tenuous interstellar plasma is a ran- 
dom medium where inhomogeneities of the density of free elec- 
trons produce fluctuations of the index of refraction. The radio 
wavefronts emitted by a point-like source (pulsar) and propagat- 
ing through such a medium produce rapid intensity fluctuations 
at the Earth. These fluctuations are due to diffraction speckles 
in the plane of the observer with a short time scale of several 
minutes and over a narrow bandwidth (~ 1 MHz) at radio fre- 
quencies. This is the DISS phenomenon (Diffractive interstel- 
lar medium scintillation). This random medium produces also 
slow intensity variations of 10-100 % over several days or more 
and over a broad bandwidth (>100 MHz). This is the RISS phe- 
nomenon (Refractive interstellar medium scintillation). Rickett 
(1988) proposes a concise presentation of these phenomena and 
of their main observables. Rickett (1977), Blandford & Narayan 
(1985) and Rickett (1990) provide a complete derivation of these 
observables as functions of the random medium physical prop- 
erties. Diffractive scintillation in dynamic spectra of pulsars is 
caused by medium inhomogeneities, whose scale so is called 
diffractive scale or coherence length. From observations of pul- 
sars at frequency around 1 GHz, this diffractive scale is between 
10 5 m and 10 8 m. Refractive scintillation of pulsars depends on 

1 Nancay radiotelescope is an instrument of the Observatoire de 
Nancay, France 



medium inhomogeneities whose scales are about the scattering 
disk radius r s = 8 S L where L is the distance to the equivalent 
turbulent screen and 6s = is the scattering angle. From 

pulsar observations, this refractive scale r$ is between 10 10 m 
and 10 13 m. These two main scales of the ionized interstellar 
medium are suggestively sketched in Cordes, Pidwervetsky & 
Lovelace (1986). The line of sight to a pulsar moves through 
the medium with a significant transverse velocity V± and con- 
tinuously samples inhomogeneities of sizes so and r$ so that the 
diffractive timescale is td = sq/V± and the refractive timescale is 
? r = LOs /Vj_. The decorrelation bandwidth Av,/ e = l/lnrs ob- 
served in dynamic spectra of pulsars is associated with the tem- 
poral broadening of pulsars t$ ~ L6 2 J2c (Rickett 1988). The 
refractive modulation index (m r - rms of intensity / mean in- 
tensity) is an important quantity directly comparable to observa- 
tions and is theoretically m r = 1.08(Ayrf/y)^ 6 (Appendix B in 
Gupta, Rickett & Coles 1993). Pulsar dynamic spectra provide 
measurements of the diffractive parameters tj and Avj c that yield 
the turbulence strength C 2 (e.g. Gupta, Rickett & Lyne 1994). 
Long term pulsar flux density series yields the refractive param- 
eter t r that is typically in tens of days and the modulation index 
m r < 1 (Stinebring et al. 2000). 

Armstrong, Rickett, & Spangler (1995) analyze a large cor- 
pus of radioastronomy observations to establish that the 3-D 
space power spectrum Pt,n(Mx, Qy, <lz) of the free electron density 
5n e (x,y,z) in the local interstellar medium (< 1 kpc) is of the 
Kolmogorov turbulence type. They give evidence of the power 

spectrum Pw(<lx> %, 9z) - C\ X ^ ^q\ + q 2 y + q 2 j valid over 

6 decades, 10~ 12 rrT 1 < Jq\ + q 2 + q\ < 10~ 6 rrT 1 . They sug- 
gest that it might extend over 12 decades. The derived turbulence 
strength by the same authors is C 2 ~ 10~ 3 irT 20 / 3 . This value 
can be considered as typical but large deviations are found de- 
pending on the celestial direction 10~ 4 < C 2 < 10 _1 m -20 ' 3 
(Johnston, Nicastro & Koribalski 1998; Cordes, Weisberg & 
Boriakoff 1985). 

The duration of ESE's (weeks to years; Maitia, Lestrade & 
Cognard 2003) and the typical transverse velocity of pulsar (tens 
of km/s; Gupta, Rickett & Lyne 1994) would make these events 
associated with large-scale inhomogeneities of free electrons, 
i.e. the refractive scales rs (10 10 m and 10 13 m) following our 
conjecture. To test this concept, we simulate the scintillation of 
the pulsar B 1937 +21 at 1.41 and 1.7 GHz with a Kolmogorov 
phase screen as large as 10 12 mx 10 12 m sampled by 131kx 131k 
pixels. The 3-D nature of this simulation is reduced advanta- 
geously to a 2-D problem in the thin screen approximation valid 
in the interplanetary and interstellar tenuous plasma. Following 
Salpeter (1967) and Lovelace (1970), this approximation is the 
limit where the ray lateral deflection is weak relative to the radial 
phase fluctuations. In this condition, the phase of a wavefront 
reduces to <f>(x, y) — J Q Ar e 6n e (x, y, z)dz, with the electron radius 

r e = 2.8179xl0~ 13 cm. The propagation of a wavefront from a 
source to an observer can be computed by the Kirchhoff-Fresnel 
diffraction integral (Born & Wolf 1999). This integral relates the 
electrical field E(x') in the observer's plane to the field E s emit- 
ted by a distant point-like source. Classically, this integral yields 
the diffraction pattern of an unperturbed aperture. We modify 
this integral by adding the Kolmogorov phase field <Pk(^) of a 
turbulent equivalent screen to model the wavefront corrugations 
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Fig. 1. Sketch of our computation of the dynamic spectrum of 
a pulsar seen through a 2-D turbulent equivalent phase screen. 
The geometric path I in the Kirchhoff-Fresnel integral, eq. (1), 
is from the point source pulsar P to the radiotelescope position 
T, by crossing the phase screen at the running pixel S (m, n). 



caused by the ionized ISM. The electric field in the observer's 
plane with the Kirchhoff-Fresnel integral is : 

E(x') = -^ fe^m^dS (1) 
^|L| Js 

where A is the wavelength, x' is the observer position vector in a 
plane parallel to the sky plane, x is the position vector of a point 
in the phase screen plane S also parallel to the sky plane, L is the 
vector between the observer and the screen planes. In addition, 
|x'| and |x| are very small with respect to |L| in our application 
below and this justifies the factor is constant in eq. (1). \(\ is 
the full geometric path from the pulsar to the observer through 
the screen. The geometry is sketched in Fig. 1. 

The 2-D power spectrum of this phase screen Pi${q x , q y ) is 
related to the 3-D Kolmogorov power spectrum of the interstel- 
lar electronic density Vw(q x , q y , q z = 0) as shown in Lovelace 
(1970) and Lovelace et al. (1970) in the thin screen approxima- 
tion : 

PiM*' = 2nz(Ar e ) 2 P 3N (q x , q y , q, = 0) (2) 

Z is the screen thickness. The thin screen approximation refers to 
a small deviation in the propagation direction of the wave rather 
than to a small thickness (Salpeter, 1967). 



In Appendix A, we derive the expression of the phase field 
<Pk(x, y) by using this relationship and the definition of the 2-D 
average power spectrum Pi^{q x , q y ) over the square screen sur- 
face N 2 Ar 2 recalled here as: 

_ , . \r 2 4qx,q y )\ 2 

P^(q x ,qy) = —r— 2 — (3) 

N 2 Ar 

T~2(f,(q x , q y ) is the Fourier Transform of the phase field 4>k(x, y). In 
Appendix A, (p K (x,y) is cast into a discrete form for computation 
by FFT The observer's plane x = (x, y) becomes a grid (m, n) 
such that x — m Ar and y = n Ar with the grid step Ar. The 
spatial frequency plane (q x , q y ) becomes the grid (k, 1) such that 
q x — k dq x and q y - I dq y . The frequency step dq x is chosen 
to be the lowest spatial frequency q m i n , x = 1 /NAr along x and, 
similarly, dq y = 1 /NAr along y. The grid step Ar is chosen to 
be so/4 after the tests described in Appendix B and adopting 
the notation so for the coherence length. The maximum spatial 
frequency q max = in the integral A. 5 satisfies the Nyquist 
sampling rate. The discrete formula reduced to the case of the 
square screen is: 

4> K (m,n) = 2n^ b Ar e JzC^ (NAr)- l+m 

k=+N/2 l=+N/2 

X J] Y J {k 2 + l 2 y m e-^ ikm+ln) e-^ l) (4) 

k=-N/2 l=-N/2 

The grid step Ar conveniently appears as part of the mul- 
tiplying factor in eq. (4) so that the double summation of this 
equation can be computed once and stored in a file for use with 
different values of Ar. This was particularly useful for the tests 
below made with several phase screens constructed by adjusting 
the multiplying factor. 

We used a complex Hermitian symmetric spectrum to make 
the phase field ^(m, n) real, i.e. The complex coefficients of 
eq. (4) are conjugate by applying tf/(k, I) = -if/(-k, -I) in the 
construction of the screen. Eq (4) is suitable for computation by 
FFT where the elements of the input complex array are the coef- 
ficients c(k, I) = -c(-k, -I) = [k 2 + fiy m e ±bKkJ >. We generate a 
random phase field <f>ic(m, n) by making the Fourier phase tf/(k, I) 
a random variable uniformly distributed over [0, 27r] (Rice 1944 
p. 287). We had to make the adjustment factor b = 2 in eq.(4) so 
that the phase structure function D^ir) computed directly from 
the screen yields a coherence length that matches the theoretical 
value sq. Coles et al. (1995) have devised a method to randomize 
a phase screen. Instead of our complex numbers be~" l ' ( - k,l) in eq. 
(4), with the random phase tff(k, I) described above, Coles et al. 
(1995) set complex random numbers made of independent zero- 
mean Gaussian random variables x and y with variance a 2 for 
their real and imaginary parts. In polar coordinates, these vari- 
ables x and y have magnitudes following a Rayleigh probability 
distribution and phases distributed uniformly over [0, 27r] as for 
our I) (see Thompson et al, 1986, p. 259-260 and reference 
to Papoulis, 1956). We have verified in generating independent 
zero-mean Gaussian random variables x and y of various cr that 
the mean amplitude of the corresponding Rayleigh probability 
distribution for cr - 1.6 is 2, i.e. the value of our b factor. Thus, 
the method we use to generate the random complex numbers in 
eq. (4) is a satisfactory approximation of the formal synthesis of 
the random phase field <p K in Coles et al. (1995). 

The intensity of the pulsar i(x', A) = |E| 2 (x', A) at wavelength 
A is computed with the discrete form of the Kirchhoff-Fresnel 
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integral where the gridded phase screen <pK(m, n) substitutes the 
continuous Kolmogorov phase 0^(x). The electric field (eq. 1) at 
the grid position p of the radiotelescope and wavelength A is : 



s structure fct c 



E(p,A) = 



AL 



2 n=H m=H 



,i[ft(p,m,n)+<l> K (m,n)] 



(5) 



where the geometric path ({p, m, n) is from the point source pul- 
sar P through the phase screen at the running pixel S (m, n) to 
the radiotelescope position T along the x-axis (p) in Fig. 1. We 
have not approximated this path by a power series in our code 
so that we assumed neither the Fraunhofer approximation nor 
the Fresnel approximation but carried out the full computation 
of €{p, m, n). The parameter dx s is the spatial resolution of the 
phase screen which does not have to be equal to the grid step Ar 
controlling the limits q min and q max of the Kolmogorov spectrum 
Pin- Since eq (5) is the discrete form of eq (1), the parameter H 
of eq (5) controls the size of the integration surface noted S in 
eq(l);S = (Hxdx s ) 2 . 

3. Construction of a Large Square Kolmogorov 
Phase Screen 

We have computed a relatively large Kolmogorov phase screen 
for the pulsar PSR1937+21 with eq. (4) on an alpha server as- 
suming turbulence strength C 2 = 10 _3 m~ 20 ^ 3 and screen thick- 
ness z to be the pulsar distance (3.6 kpc) (this assumption is dis- 
cussed below). The size of the screen is N = 2 17 , i.e. there are 
~ 131k x 13 lk pixels. Obviously, this computation could not be 
done in a single 2-D FFT and we had to resort to multiple 1-D 
FFT's. Each line of the input complex array, coefficients c(k, I) 
of eq.(4), was first transformed by a 2 17 point 1-D FFT along the 
x-axis. Then, each resulting column was transformed by a 2 17 
point 1-D FFT along the y-axis. The construction of this large 
square phase screen took 10 days on our workstation and needs 
68 GigaBytes of disk space to store all the phases but the peak 
storage during the computation was 200 GigaBytes. 

With the assumed turbulence strength C 2 = 10~ 3 irT 20/3 and 
thickness z = 3.6 kpc, the coherence length in the screen is 
so = 2.66 1 7 m at 1.41 GHz as calculated from the phase 
structure function D^sq) = 1 Rd 2 with the theoretical expres- 
sion D^is) = %nr 2 A 2 C^zf{a){a + l)~ l s", where numerically 
f(a) = 1.12 for a-5/3 (Armstrong, Rickett & Spangler 1995). 
To qualify our phase screen, we have computed the phase struc- 
ture functions D^(s) along multiple x-lines through the screen 
and, orthogonally, along y. From the log-log plots of these func- 
tions (examples in Fig. 2), we made linear least-square-fits and 
found slopes of 1.68 ±0.10 for x and 1.64±0.1 fory. The theoret- 
ical slope of D^s) is a — 5/3 for an isotropic turbulent medium 
with the exponent (3 = 1 1/3 in P 3N . We have also measured the 
coherence lengths from these phase structure functions: r con = 
2.4 + 0.6 10 7 m along x and r con = 3.2 + 0.6 10 7 m along y. These 
values are consistent with the value so = 2.66 1 7 m predicted 
by the structure function D^(so) = 1 Rd. This is also consistent 
with the value obtained from the equation r con = V/td -2.1 10 7 
m, where V — 50 km/s and ta = 7 min is the observed diffractive 
timescale toward the pulsar B1937+21 (Ryba 1991). We note 
that it is coincidental that so has the expected value by using the 
conventional value of C 2 for PSR1937+21 and our assumption 
above that the thickness z equals the pulsar distance. 

4. Results: Simulation of the Scintillation of the 
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Fig. 2. Examples of the phase structure functions D^x) and 
D$(y) computed for 1-D cuts along the x and y axis across 
the whole 131k x 131k phase screen <pK{m,ri) constructed with 
C 2 = 10~ 3 , z — 3.6 kpc and Ar = so. The slopes of these struc- 
ture functions measured from the log-log plots are satisfacto- 
rily close to the theoretical value a = 5/3 for /3 = 11/3 of the 
Kolmogorov spectrum P 3N . 



Pulsar B1 937+21 at 1.41 and 1.7 GHz 

Prior to the full scale simulations, we studied the convergence 
of the computation of the dynamic spectra that depends on four 
parameters : the size S of the integration surface of eq. (1), the 
spatial resolution dx s to read the phase screen file, the grid step 
Ar in the phase screen and the size N of the phase screen (eq. 4) 
which control the minimum and maximum spatial frequencies 
qmin = 77X7 and q max = f the Kolmogorov spectrum P 3N . 
The parameter N was fixed to 2 17 , limited only by the maximum 
disc space available with our computer to store the phase screen, 
and this corresponds to a period of 6 months if the phase screen 
drifts across the observer line of sight at V = 50 km/s. We have 
empirically determined the other three parameters by a sequence 
of tests that are described in detail in Appendix B. 

In Fig. 3, we present the intensity of the pulsar B1937+21 
simulated during this period of six months at 1.41 and 1.7 GHz. 
The intervening phase screen (pnim, n) is constructed with C 2 = 
10" 3 m _20/3 and z = 3.6 kpc that are the parameters of B1937+21 
(Ryba 1991; Johnston, Nicastro & Koribalski 1998), N = 2 17 
and Ar = so/4 in eq. (4) based on our conclusion in the 
second test (Appendix B). The coherence length in the phase 
screen is 2.6 10 7 m thus the lower and upper limit scales of the 
Kolmogorov spectrum are q„ 

and q„ 



imin - 217x2.6 10 7 ~ 3.4 10 12 

2x2.6 io? m ■ ^ e adopt the integration surface 
S = (4rj) 2 to compute eq. (5) from the conclusion in the first 
test (Table 1). We use the resolution dx s = so/4 (Table 2) to read 
the phase screen file into eq. (5). This resolution is conservative 
relative to our conclusion in the third test; it was dictated for 
algorithmic peculiarities to the expense of computing speed. A 
dynamic spectrum was computed in these conditions every other 
1.25 days. Averaging these dynamic spectra provides the inten- 
sity measured by the telescope over the integration time 70 min 
and the bandwidth 8.8 MHz, sampled over 32 x 32 time bins 
and frequency channels in our computation. The velocity of the 
screen is V = 50 km/s. This simulation took twice 1 .5 months 
on our alpha server. 
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Fig. 3. Intensity of B 1937+21 simulated at 1.41 GHz and 1.7 GHz. The intervening phase screen is constructed with C 2 , = 10 -3 , 
z = 3.6 kpc, N = 2 17 and Ar = sq/4 in eq. (4) and S = (4r s ) 2 and dx s = sq/4 in eq. (5). Between days 30 and 40, there is an event 
occurring simultaneously at 1.41 and 1.7 GHz which is alike the behavior of the flux density recorded in direction of B1937+21 
at Nancay in October 1989 and interpreted as an Extreme Scattering Event. The mean intensities < / > are 1.46 and 1.60 and the 
modulation indexes m are 0.49 and 0.67 for these two series. 



Overall, Fig. 3 shows the typical behavior of flux density 
recorded for pulsars; for instance B1937+21 at Nangay at 1.41 
and 1.7 GHz in Figs. 1 and 2 of LRC98. In addition, there is the 
interesting feature labeled "ESE" and present simultaneously at 
1.41 and 1 .7 GHz in Fig. 3 from day 30 to day 40. The significant 
drop of flux density (60 % below the mean value) and the low 
rms during this ESE period (rms=0.10 during the ESE versus 
rms=0.71 off ESE) for these 10 days is alike the ESE observed 
toward B 1937+21 in October 1989 (Cognard et al. 1993). In this 
previous paper, we suggested this phenomenon was caused by 
a discrete cloud of plasma following the standard interpretation 
of Fiedler at al. (1987). Our simulation shows instead that this 
event can arise naturally because of the turbulence in the ionized 
interstellar medium as conjectured. We have labeled this event 
"ESE" in Fig. 3. The ratio of the duration of this event to the 
length of our simulation is ~ 10days/180days ^ 5%. Although 
this might be fortuitous, this percentage is consistent with the 
observations of B 1937+21 at Nancay for which the rate of oc- 
currence of "ESE's" is 4% (LRC98). Figure 4 shows the cross 
correlation function of the flux densities at 1.41 and 1.7 GHz 
(Figure 3). The zero-lag value ^0.8 is consistent with the ob- 
served value 0.93+0.05 (LRC98). 

The duration of the event labeled "ESE" in our simulation 
is comparable to the event observed at Nancay in October 1989 
(12 days) but is short compared to three others of duration 1-3 
months identified in the radio light curves of B1937+21 simul- 



Lag (day) 



Fig. 4. Cross-correlation function of the flux densities at 1.41 
GHz and 1.7 GHz of Figure 3. 



taneously at 1.41 and 1.7 GHz by LRC98. The extraordinary fo- 
cusing of pulsar waves is caused by the large scale fluctuations 
of the electronic density in the interstellar medium. The largest 
fluctuation scale in our simulation is 2 17 x 2.6 10 7 /4 ~ 10 12 m 
while the ionized interstellar medium is made of fluctuations up 
to 10 18 m (Armstrong, Rickett and Spangler 1995). We expect 
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that longer "ESE"s require to include fluctuations of larger scales 
in our simulation. We plan to extend this calculation on a super- 
computer. We note here that the relatively long events labeled 
"ESE" in Hamidouche, Lestrade, and Cognard (2002) are from 
a simulation that was improperly done with an integration sur- 
face as small as S = (0.7 r s ) 2 chosen before we made the tests of 
the present paper. 

Finally, we note that the anti-correlation coefficient -0.62 be- 
tween flux density and TOA in the simulation conducted with 
a rectangular phase screen is reported in Hamidouche (2003), 
and is comparable to the correlation coefficient -0.43 to -0.73 de- 
pending on data segments of the Nancay observations in LRC98. 
We have not calculated this correlation in the final simulation. 

5. Conclusion 

We have presented the simulation of the scintillation of pulsars 
carried out within the frame of Physical Optics by extending the 
seminal work by Cordes, Pidwervetsky & Lovelace (1986). We 
have built a large square Kolmogorov phase screen of 131k x 
131k pixels on an alpha server and shown how to set properly 
the parameters of the Kirchhoff-Fresnel integral by conducting 
several tests. Our final simulation was done for the condition 
of turbulence known in the direction of the pulsar B1937+21 
(C 2 ~ 10~ 3 and distance z — 3.6 kpc) and the radio light curves 
were generated at 1.41 and 1.7 GHz for a period of 6 months. 

These two light curves exhibit simultaneous variations at 
1.41 and 1.7 GHz that are alike the "ESE" observed at these fre- 
quencies at Nancay in October 1989 and that lasted ~ 10 days. 
Our simulation shows that this observed event can be caused 
naturally by the turbulence in the ionized interstellar medium 
instead of invoking the crossing of a discrete over pressured ion- 
ized cloud on the line of sight as in the model of Fiedler et al. 
(1987). We think that longer events could occur by including the 
electronic density fluctuations of larger scales in the construction 
of the Kolmogorov phase screen used in the simulation. 

Finally, we note that Deshpande (2000) stresses that the 
opacity differences in HI and other species measured over a 
transverse separation xo result from all scale of the 3-D power 
spectrum of the opacity fluctuations while they are currently, 
and erroneously in his opinion, interpreted as over pressured 
and overdensed cloudlets of size xo- Our simulation of "ESE's" 
caused by the turbulent ionized interstellar medium strengthens 
this opinion for the neutral gas. 
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Appendix 



Appendix A: Construction of a Kolmogorov Phase 
Screen 

The definition of the 2-D phase power spectrum of the phase 
field (f>(x, y) averaged over the rectangular screen surface NMAr 2 
is : 

_ , . \T2*(qx,q > )\ 2 

Pi4,(qx, q y ) = 2 (A. 1) 

NMAr 

where T2<f,(q x , qy) is the Fourier Transform of <f>(x,y). Hence, the 
module of the Fourier component for the frequencies (q x , q y ) is: 

\r^(q x ,q y )\= ^NMAr 2 P 2lP (q x ,q y ) (A.2) 
Lovelace (1970) and Lovelace et al. (1970) have established: 

P2</,(q x , q y ) = 2nz(Ar e f¥ m (q x , q y , q z = 0) (A.3) 

Hence, the complex Fourier components of the phase field 
<p(x, y) are: 

r^(q x ,q y ) = -fhtXr e ^NMAr 2 zC 2 (q 2 x + q 2 y Y^ 2 e** 9 >«> ) (A.4) 

that must satisfy the complex Hermitian symmetric spectrum 
with the Fourier phase i/f(q x , q y ) = -tf/(-q x , -q y ) to make the 
phase field (p(x, y) real. tff(q x , q y ) is a random variable uniformly 
distributed over [0, 2n] as prescribed by Rice (1944, p. 287). C 2 
is the turbulence strength and z is the propagation length. 
The phase field (f>(x, y) can be computed by the inverse Fourier 
Transform of T74 (q x , q y ). Although is a power-law, the 
Fourier integral is finite because there is an outer scale in the 
ionized interstellar medium that makes P^ N becomes zero for 
small q rather than infinity when q — > 0. Also, we point out that 
there is the factor {2nY 2 in this integral to be consistent with 
the Fourier Transform definition used to define in Rickett 
(1977, eq. 6) and used to demonstrate the Lovelace relationship 
(Lovelace 1970, eq. 36). This inverse Fourier Transform is: 

<p{x,y) = (2ny 2 ^-Ar e Ar^NMzC 2 N x 

(ql + q]r PIA e- mq '* } e- i{q > x+q ^dq x dq y 



This integral was split into two parts to avoid the singularity of 

the power-law function. In eq. (A.5) for the rectangular screen 

case, the minimum frequencies are q x , mm = j^p and q yM i„ 
i 
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Appendix B: Tests for the Computation of the 
Dynamic Spectra 

Prior to the full scale simulations, we studied the convergence 
of the computation of the dynamic spectra that depends on four 
parameters : the size S of the integration surface S of eq. (1), 
the spatial resolution dx s to read the phase screen file, the grid 
step Ar in the screen phase, and the size N of the phase screen 
(eq. 4). These latter two parameters control the minimum and 
maximum spatial frequencies q min = ^ and q max = ^ of the 

Kolmogorov spectrum Pw- The parameter N is fixed to 2 17 in 
our computation and we have empirically determined the other 
three parameters. 

As a first test, we study the effect of S of eq. (1) while keep- 
ing the parameters Ar and dx s fixed and equal to so. We express 
the dimension of the integration surface S in terms of the scat- 
tering radius r s since, as shown below, the required S to make 
the dynamic spectra numerically convergent is a few times this 
scale. In fact, Coles et al. (1995) have demonstrated that the in- 
tegration surface in the phase screen has to be larger than the 
scattering disk. The physical idea is that the scattering disk is 
the approximate area in the phase screen irradiating a single lo- 
cation in the observing plane. We chose the following sizes S 
for our test: S = (lr s ) 2 , S = (2r s ) 2 , S = (4r s ) 2 , S = (6r s ) 2 , 
S = (8r s ) 2 , S = (10r s ) 2 . For each size S, we compute the whole 
series of dynamic spectra at 1.41 GHz by shifting S by 2.5 days 
across the screen using V=50 km/s. Each dynamic spectrum is 
70 min x 8.8 MHz in size, similarly to observations of B 1937+21 
at Nancay, and sampled over 8x8 pixels. For each series, we de- 
rive the mean intensity < I > and the modulation index m = ^ 
in Table B.l. This table shows that these two indicators converge 
when S reaches (4rj) 2 . In Fig. B.l, we show only three series 
for clarity; S = (2r s ) 2 , S = (4r s ) 2 and S = (10r s ) 2 to depict the 
convergence process. In addition, we illustrate this convergence 
process in Fig. B.2 by showing dynamic spectra computed at 
the same position in the observer plane but for the 6 integra- 
tion surfaces mentioned above. From this first test summarized 
by Table B.l, Figs. B.l and B.2, we conclude that convergence 
of dynamic spectra starts for S = (4r s ) 2 since improvement in 
the average intensity and index m are less than 10 % for larger 
surfaces S . 



Table B.l. Sensitivity of the averaged intensity and modulation 
index to the integration surface at 1 .4GHz a 



Integration 


Nber of 


< / > 


Index 1 


surface 5 


intensities 






(ir s ) 2 


303 


0.28 


0.95 


(2r s ) 2 


302 


0.83 


0.80 


(4r s ) 2 


299 


1.39 


0.63 


(6r s ) 2 


296 


1.52 


0.59 


(8r,) 2 


292 


1.55 


0.58 


(Wr s ) 2 


289 


1.56 


0.58 



Ar = dx s = s kept fixed. 



As a second test, we study the impact of the grid step Ar in 
the screen phase while keeping fixed S = (4rj) 2 and dx s = 
so. The parameter Ar is directly related to the minimum and 



maximum spatial frequencies q min = and q max = ^ of 
the Kolmogorov spectrum Pj, N as already mentioned. Note that 
Ar appears conveniently in the multiplying factor of the phase 
screen (picim, n) in eq. (4) and hence it can be easily changed to 
any value to modify q max and, correlative ly, q min . The parame- 
ter N is not amendable after the double summation of eq. (4) 
has been calculated and stored in a computer file. In principle, 
we would have liked to tune N in order to keep q min = ^ 
unchanged while Ar is adjusted to increase q max = One 
expects that the high frequency part of the Kolmogorov spec- 
trum becomes insignificant in shaping the dynamic spectrum 
when q max is sufficiently high while the low frequency part mod- 
ifies it drastically. This test is difficult to implement in practice 
because it requires computing several screens with different N 
while keeping the same random phases if/{k, I) for the lower part 
of the spectrum. Instead of this approach, we simulate the effect 
by superimposing a corrugated surface of period q^ 1 and am- 
plitude §<p upon the original screen. The aim is to seek which 
perturbating surface (q, 8(p) degrades significantly the dynamic 
spectrum when compared to the one computed with the original 
phase screen. This test covers the following corrugated surfaces 
: 8<f> = 101° and q = 5^ ; 6<p = 57°and q = 5<f> = 32° 

and 9 = 2^72' 6 <t> = 18 ° and 9 = 5^74= 6< t> = 10 ° for 9 = 2^78- 
These amplitudes §<p are derived from the phase structure func- 
tion D^s) for the separations s = 2 sq, s = so, s = so/2, 
s = so/4 and s = so/8. Fig. B.3 shows the comparison of the dy- 
namic spectra computed with these five perturbations (q, d<f>) as 
well as the dynamic spectrum simulated with the original screen 
(dip = 0) constructed with C 2 = 10~ 3 , z = 3.6 kpc, N = 2 17 and 
Ar = so/4. The morphology and averaged intensities of these 
dynamic spectra indicate that perturbations are significant above 
(8<p = 32°, q = 2^/2X i- e - qmax must be at least as high as 5^72* 
i.e. Ar < so/2, for convergence of the computation. 
We complement this test done at a single position of the observer 
by simulating three intensity series across the full screen with the 
grid step Ar = dx s set to lso, so/2 and so/4. The dynamic spectra 
are sampled over 8 x 8, 16 x 16, 32 x 32 pixels, respectively, to 
synthesize the same 70 min x 8.8 MHz domain. Table B.2 shows 
that the mean intensity < / > and modulation index m are within 
10 % over these three cases. From this second test we conclude, 
somewhat conservatively, that Ar = so/4 is necessary for the 
convergence of the computation, i.e. the Kolmogorov spectrum 
P 3N must include q max = jj^. 



Table B.2. Sensitivity of the averaged intensity and modulation 
index to q max a 



Ar 


Qmax 2 Ar 


Sampling 


Nbof 


< / > 


Index m 






of dyn. sp. 


intensities 






so 


1 

2 so 


8x8 


299 


1.39 


0.63 


s /2 




16x 16 


142 


1.43 


0.54 


so/4 


2 sp/2 
2~W4 


32x32 


69 


1.58 


0.50 



a Ar = dx s = s kept fixed. 



As a third test, we study the impact of the spatial resolution dx s 
used to read the phase screen file into eq. (5). Although we have 
just concluded that Ar = s /4 is necessary to include enough 
high frequencies, the phases can possibly be read with a lesser 
resolution. Using the phase screen file constructed with C 2 = 



M. Hamidouche and J.-F. Lestrade: Simulation of the interstellar scintillation and the extreme scattering events of pulsars 



9 



1(T 3 , z = 3.6 kpc, N = 2" and Ar = s /4, we test the three 
cases: dx s = sq in reading every other four phases, dx s = sq/2 
in reading every other two phases and dx s = *o/4 in reading 
every phase. The integration surface is fixed to S = (4r^) 2 in 
this test. Fig. B.4 shows the three dynamic spectra corresponding 
to these resolutions. As the Figure shows, they have the same 
morphology and their averaged intensities changed by only 6 %. 
We extended the resolution further to dx s = 2sq and found then 
that the morphology of the dynamic spectrum becomes abruptly 
very dissimilar (not shown in Fig. B.4) meaning it had not yet 
converged. We also computed the intensity series for significant 
fractions of the screen and report < / > and m in Table B.3. From 
this third test, we conclude that dx s = so is sufficient. 

Table B.3. Sensitivity of the averaged intensity and modulation 
index to dx s a 



Resolution dx s 


Sampling 


Nbof 


<I> 


Index m 




of dyn. sp. 


intensities 






Ls'o 


8x8 


69 


1.57 


0.54 


s /2 


16 x 16 


31 


1.87 


0.49 


so/4 


32x32 


18 


1.98 


0.49 



S = (4r s ) 2 and q max = ^ kept fixed. 



Finally, we note that the resulting modulation index m (~ 0.5) 
is larger than the theoretical value m r = 0.32 derived for 
purely refractive scintillation by Gupta, Rickett & Coles 1993 
(m r = 1.08 (^-) 1/6 ). This is because the size of the dynamic 
spectra simulated (70 min x 8.8 MHz) is not much larger than 
the diffractive patches f^xAvrfc = 8.8minxl.l MHz (td = sq/V 
and Av dc = l/(2nT S )). This size (70 min x 8.8 MHz) we used 
to average out diffractive scintillation is not large enough and so 
we have a blend of both refraction and diffraction into m. We 
note also that the normalized intensity < / > is ~ 1 .5 rather than 
unity but this is a statistical fluctuation. For instance, we found 
< / >= 0.77 and m = 0.51 for scintillation simulated along a 
full Ar s wide track in another part of our large phase screen. In 
Fig. B.5, we provide an example of 1 1 dynamic spectra 2.5 days 
apart and their summed autocorrelation function. The visibility 
at o- = s (~ 60%) in < E(s)E*(s + <x) >= exp[-0.5((r/ s Q ) a ] 
(Rickett 1988) is closely delineated by the contours 65% in the 
2-D autorrelation function of Fig. B.5. The half widths of this 
function are ~ 4 min and ~ 0.5 MHz and are consistent with the 
theoretical values td = 8.8 min and Ava c = 1.1MHz. 
In Fig. B.6, we compute the normalized autocorrelation func- 
tion (flc/(r)/flc/(0)) of the time series of the intensity shown 
in Fig. B.l for the cases S = (4-r s ) 2 . This function shows the 
slow decline, ~ 10 days at half-maximum of acf(l)/acf(Q), ex- 
pected because of the correlation induced by the long refrac- 
tive scales of the phase screen. This is consistent with the the- 
oretical refractive time scale t r = L9 S /V (Rickett 1988) of 15 
days with L = 1.8 kpc, V = 50 km/s and the scattering angle 
6 S = = 1.25 10~ 9 Rd, with A = 0.21 m and s = 2.66 10 7 m. 
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Fig. B.l. Test on the parameter S : simulations of the intensity of a pulsar at 1.41 GHz every other 2.5 days observed through the 
phase screen <p^(m,ri) constructed with C 2 , = 10~ 3 , z = 3.6 kpc, N = 2 17 and Ar = sq. The spatial resolution used to read the 
screen file is dx s = Ar. The three light curves are for the integration surfaces S = (2rs) 2 (filled circle), S = (4rs) 2 (filled triangle), 
S = (10r s ) 2 (empty circle). This test demonstrates that the computation of the dynamic spectra has converged for the integration 
surface S = (4r s ) 2 . 
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Dynamical spectrum S = ( 1 r s ) z 
averaged intensity = 0.37 



Dynamical spectrum S = ( 2r g ) z 
averaged intensity = 2.28 




Dynamical spectrum S = ( 4r s ) z 
averaged intensity = 2.44 




Dynamical spectrum S = ( 6r s ) z 
averaged intensity = 2.57 




Dynamical spectrum S = ( 8r s ) 2 
averaged intensity = 2.59 




Dynamical spectrum S = (10r s ) 2 
averaged intensity = 2.61 




Fig. B.2. Test on parameter S : dynamic spectra computed, at 1.41 GHz, at the same telescope location but for different sizes S of 
the integration surface S in eq. (1). The six cases shown are for S = (lr s ) 2 , (2r s ) 2 , (4r s ) 2 , (6r s ) 2 , (8r s ) 2 , (10r s ) 2 . The total size 
of the dynamic spectrum is 70 min x 8.8 MHz sampled over 8x8 pixels, i.e. unity corresponds to 1.1 MHz along the frequency 
axis (vertical) and to 8.8 minute along the time axis (horizontal) with the screen speed V = 50km/s. The phase screen (p K (m,n) is 
constructed as in Fig. 2. This Figure shows that the dynamic spectrum converges from the surface size S = (4r s ) 2 . 




Dynamical spectrum i 
averaged intensity = 



= 32° 
1 .6792", 




Dynamical spectrum (50 = 57° ) 
averaged intensity = 0.94833 




Dynamical spectrum (60=101°) 
averaged intensity = 0.24807 




Fig. B.3. Test on parameter Ar: dynamic spectra computed at 1 .41 GHz at the same telescope location and with the same integration 
surface size S = (4r s ) 2 but with different grid step Ar for the phase screen 0#(m, n) in eq. (4). The limits of the Kolmogorov spectrum 
are directly related to Ar and N as q m i„ = ^ and q max = ^ • A full test would require to construct effectively several Kolmogorov 
phase screens in decreasing Ar to modify q max and in increasing N to keep the same q min . This is very difficult to implement 
numerically (see text) and, instead, we have simulated this effect in superimposing a corrugated surface onto the original phase 
screen. This surface is characterized by an amplitude &<p for some spatial frequency q. The idea is to find out which perturbating 
surface 8(p is required to modify significantly either the averaged intensity or the morphology of the resulting dynamic spectrum. 
Our computation covered the following cases : 5(p = 10° for q = jj^', 54> = 18° for q = 5^74 ^ = 32° for q = 2 ^ /2 ; 6(f> = 57° 
for q = j-L; 5(f> = 101° for q = j^. These amplitudes 



s = l/8so, s = 1/4 so, s 



4 J ' 

1/2j , s = l*o, * : 



are derived from the phase structure function D^is) for the separations 
2sq, respectively. The first panel of the Figure shows the dynamic spectrum computed 



with the original Kolmogorov screen (6<f> = 0°) constructed with C 2 n = 10 3 , z = 3.6 kpc, N = 2 17 and Ar = sq/4. This test shows 



that q„ 



2 Jo /2 



i.e. Ar = sq/2, is necessary to include all relevant Fourier components for the computation of dynamic spectra. 



The gray scale and contours are the same for all panels. Unity corresponds to 1 . 1 MHz along the frequency axis (vertical) and to 8.8 
minutes along the time axis (horizontal) with the screen speed V = 50km/s. 
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Dynamical spectrum (dx =s q /1) Dynamical spectrum (dx s = s Q /2) Dynamical spectrum (dx s = s /4) 



averaged intensity = 1.68 averaged intensity = 1.75 averaged intensity = 1.78 




2 4 6 8 5 10 15 10 20 30 



Fig. B.4. Test on the parameter dx s : dynamic spectrum computed at 1.41 GHz at the same telescope location and with the same 
integration surface S = (4rs) 2 but with different spatial resolution dx s in reading the phase screen file. The phase screen <pK{tn, n) 
has been constructed with C\ = 10~ 3 , z = 3.6 kpc, N = 2 17 and Ar = so/4. Only every other 4 pixels are read from the phase file 
into eq. (5) in the case dx s = so; only every other 2 pixels for dx s = So/2; every pixel for dx s = sq/4 matching Ar = sq/4 in that 
case. The size of the dynamic spectrum is the same for the three panels, 70 min x 8.8 MHz, sampled over 8 x 8, 16 x 16 and 32 
X 32 pixels, respectively, and unity is correspondingly 1.1 MHz, 0.55 MHz, 0.275 MHz along the frequency axis (vertical) and 8.8 
min, 4.4 min, 2.2 min along the time axis (horizontal) with the screen speed V = 50km/s. This test shows that the spatial resolution 
dx s = so is sufficient to make the computation of dynamic spectra and averaged intensities < i(x', A) > convergent. 
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Fig. B.5. Sample of simulated dynamic spectra at 1.4GHz and their summed autocorrelation function (last panel). All the dynamic 
spectra are plotted with the same intensity scale, they are 1.25 days apart with the adopted screen speed V = 50 km/s. They show 
a high variability in morphology and averaged intensity as expected for pulsars. The structure in these dynamic spectra is similar 
to speckles seen in images at optical wavelengths. Their summed autocorrelation provides a mean to measure the diffactive time 
scale and decorrelation bandwidth. We found they match the theoretical values (see text). Unity corresponds to 0.275 MHz along 
the frequency axis (vertical) and to 2.2 minutes along the time axis (horizontal) for both the dynamic spectra and the summed 
autocorrelation function. This simulation is done with C\ = 10~ 3 , z = 3.6 kpc, dx s = Ar = sq/4, N = 2 17 and S = (4rs) 2 . 




B.6. Autocorrelation function of the time series of the intensity computed with S = (4-rs) 2 and shown as filled triangles in 
B. 1 . The refraction time scale of ~ 10 days at half-maximum of this function is consistent with the theoretical value of 15 days. 



